logo

# Libraries
library(tidyverse)
library(sf)
library(USAboundaries)
library(ggthemes)
library(cowplot)
library(sp)
library(leaflet)



Measurement threshold (1960 - 2000)

A minimum threshold of 10 distinct measurements with at least 1 measurement every 5 years was applied to filter for only the regularly monitored wells. investigate regularly monitored wells




# Data from USGS National Water Information System (NWIS)

# Arizona state shape (CRS = 5070)
az = us_states() %>% 
  filter(name == 'Arizona') %>% 
  st_transform(5070)

# add 'well #' for each unique well
usgs = az_nwis_unique_sites
usgs$wellid = paste("well", 1:nrow(usgs))

# remove wells located at same lat/long but have different well ID
usgs_spatial = usgs %>% 
  st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>% 
  st_transform(5070) %>%
  as_Spatial() %>% 
  remove.duplicates() %>% 
  st_as_sf() %>% 
  group_by(wellid)

# match 'well #' from spatial data frame to corresponding well in time series data frame
usgs_time = az_nwis_all %>% 
  group_by(site_id) 
usgs_time = left_join(usgs_time, select(usgs, site_id, wellid), by = "site_id") 
usgs_time = usgs_time %>% filter(wellid %in% usgs_spatial$wellid)

# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh = usgs_time %>% 
  group_by(wellid) %>% 
  filter(measurement_dist >= 10) %>% 
  mutate(measure_period = year - lag(year)) %>% 
  filter(any(measure_period > 5))
usgs_time = usgs_time %>% filter(!wellid %in% thresh$wellid, measurement_dist >= 10) %>% 
  mutate(measure_period = year - lag(year))
usgs_spatial = usgs_spatial %>% 
  filter(!wellid %in% thresh$wellid) %>% 
  filter(measurement_dist >= 10)


# Data from Arizona Department of Water Resources (ADWR)

# add 'well #' for each unique well
state = adwr_unique_sites %>%
  filter(date > 1960-01-01, date < as.Date.character('2000-01-01'))
state$wellid = paste("well", 1:nrow(state))

# remove wells located at same lat/long but have different well ID
state_spatial = state %>%
  st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>%
  st_transform(5070) %>%
  as_Spatial() %>%
  remove.duplicates() %>%
  st_as_sf() %>%
  group_by(site_id)

# match 'well #' from spatial data frame to corresponding well in time series data frame
state_time = adwr_all %>%
  group_by(site_id)
state_time = left_join(state_time, select(state, site_id, wellid), by = "site_id")
state_time = state_time %>% filter(wellid %in% state_spatial$wellid)

# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh_state = state_time %>%
  group_by(wellid) %>%
  filter(measurement_dist >= 10) %>%
  mutate(measure_period = year - lag(year)) %>%
  filter(any(measure_period > 5))
state_time = state_time %>% filter(!wellid %in% thresh_state$wellid, measurement_dist >= 10) %>%
  mutate(measure_period = year - lag(year))
state_spatial = state_spatial %>%
  filter(!wellid %in% thresh_state$wellid) %>%
  filter(measurement_dist >= 10)





Wells deeper than 500ft





500 - 1000ft







500 - 600ft


I decided to include these increments of 100ft as it is visually easier to track individual wells






600 - 700ft






Wells with over 50 measurements





Phoenix AMA





Well 5284





Well 5324






Harquahala AMA




Well 292






Well 504


Located just north of Harquahala AMA boundary






Kingman




Well 1139







Outside Tuscan AMA




Well 575






Wells deeper than 1000ft




Flagstaff


Well 1133





Well 1215






Well 3450


Well 3450 was monitored continuously since ~1965





Well 1242






Well 1226


Wells 1226 & 2710 had measurements recorded only for a couple years in the late 1990s






Well 2710







Between Flagstaff and Kingman




Well 5155


Northwest of Flagstaff and northeast of Kingman





Northeast of Flagstaff


Well 1669


shows a continual decline in DTW from ~1975 - 2000





Well 1411


Well 1411 is nearby to well 1669 and shows the same trend as 1669, 1411 has a max depth of 830 ft, just under 1000ft marker





Arizona Department of Water Resources (ADWR)




Using the same threshold applied to the USGS data of 10 distinct measurements with at least 1 measurement every 5 years, the total number of State monitored wells was significantly reduced.




300 - 400 ft




400 - 500 ft




500 - 600 ft




600 - 700 ft




Functions


# PLOTS A BUFFER AROUND SPECIFIED WELL AND LOCATES OTHER WELLS INSIDE THE BUFFER AREA
buffer_fun = function(df, well, buff, state) {
  buffer = st_buffer(df[well,], buff)
  near1 = st_intersection(df[,], buffer) %>% filter(measurement_dist >= 10)
  
  plot = ggplot() + 
    geom_sf(data = state) +
    geom_sf(data = buffer, fill = NA) + 
    geom_sf(data = near1, col = "red", size = .5) + 
    labs(caption = paste(nrow(near1), 'wells')) +
    theme_void() +
    theme(plot.caption = element_text(size = 22, face = "bold", hjust = 0.5))
  print(plot)
  return(near1)
}

# RETURNS DATAFRAME OF WELLS IN DTW RANGE (MIN - MAX)
dtw_range = function(df, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    mutate(sd = sd(dtw)) %>% 
    arrange(desc(date))
  return(df)
}


# PLOTS AN INDIVIDUAL WELL BY WELL ID (DTW TIME SERIES)
plot_well = function(df_time, id) {
  wells = df_time %>% filter(wellid == paste('well', id))
  plot = ggplot(data = wells, aes(x = date, y = dtw_ft)) +
    geom_line(aes(y = dtw_ft, col = wellid), size = 2) +
    scale_y_reverse() +
    labs(title = paste('Well', id),
         caption = paste('Min depth:', min(wells$dtw_ft), '\nMax depth:',
                          max(wells$dtw_ft), '\nNumber of measurements:',
                          wells$measurement_dist),
         x = 'Year',
         y = 'DTW (ft)') + 
    theme_bw() +
    theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5),
          axis.text.x = element_text(color="black", size=14), 
          axis.text.y = element_text(color="black", size=14), 
          axis.title.x = element_text(face = 'bold', color="black", size=16), 
          axis.title.y = element_text(face = 'bold', color="black", size=16),
          plot.caption = element_text(face = 'bold', color = 'black', size = 14),
          panel.grid.major = element_line(colour = "#808080"),
          panel.grid.minor = element_line(colour = "#808080", size = 1),
          legend.position = 'none') 
  print(plot)
}

# PLOTS MULTIPLE WELLS IN DTW RANGE
multi_well_plot = function(df, df_time, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    arrange(desc(date))
  df_time = df_time %>% filter(wellid %in% df$wellid)
  plot1 = ggplot(data = df_time, aes(x = date, y = dtw_ft)) +
    geom_line(aes(y = dtw_ft, col = wellid), size = 1) +
    scale_y_reverse() +
    labs(x = 'Year',
         y = 'DTW (ft)') + 
    theme_bw() +
    theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5), 
          axis.text.x = element_text(color="black", size=14), 
          axis.text.y = element_text(color="black", size=14), 
          axis.title.x = element_text(face="bold", color="black", size=16), 
          axis.title.y = element_text(face="bold", color="black", size=16), 
          panel.grid.major = element_line(colour = "#808080"),
          panel.grid.minor = element_line(colour = "#808080", size = 1),
          legend.position = 'none')
  print(plot1)
  return(df)
}


# PLOTS WELLS DTW WITHIN A SPECIFIED RANGE, COLOR GRADIENT BY DEPTH, RETURNS CORRESPONDING DATA FRAME
map_dtw = function(state, df, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    mutate(sd = sd(dtw)) %>% 
    arrange(dtw)
  plot1 = ggplot() + 
    geom_sf(data = state, size = 1, col = 'black') +
    geom_sf(data = df, aes(col = dtw), size = 1.5) + 
    theme_void() +
    theme(plot.title = element_text(face = "bold", color = "black", size = 18),
          plot.subtitle = element_text(size = 12),
          plot.caption = element_text(face = 'bold', size = 12), 
          legend.title = element_text(color="black", size=14), 
          legend.text = element_text(color="black", size=14)) +
    scale_colour_gradient() +
    labs(title = paste('DTW', min, '-', max, 'ft'),
         fill = 'Number of wells', 
         col = 'Depth to Water (ft)')
  print(plot1)
  return(df)
}